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ABSTRACT 

We study the nature of non-axisymmetric dynamical instabilities in differentially ro- 
tating stars with both linear eigenmode analysis and hydrodynamic simulations in 
Newtonian gravity. We especially investigate the following three types of instability; 
the one-armed spiral instability, the low T/|VF| bar instability, and the high T/|W| bar 
instability, where T is the rotational kinetic energy and W is the gravitational poten- 
tial energy. The nature of the dynamical instabilities is clarified by using a canonical 
angular momentum as a diagnostic. We find that the one-armed spiral and the low 
T/|W| bar instabilities occur around the corotation radius, and they grow through 
the inflow of canonical angular momentum around the corotation radius. The result 
is a clear contrast to that of a classical dynamical bar instability in high T/|W|. We 
also discuss the feature of gravitational waves generated from these three types of 
instability. 

Key words: gravitational waves - hydrodynamics - instabilities - stars: evolution - 
stars: oscillation - stars: rotation 



1 INTRODUCTION 

Stars in nature are usually rotating and may be subject 
to non-axisymmetric rotational instabilities. An analytically 
exact treatment of these instabilities in linearized theory 
exists only for inco mpressible equilibrium fluids in New- 
tonian gravity (e.g.. IChandrasekharl Il969t iTassoull Il978t 
IShapiro fc Teukolskvll983r) . For these configurations, global 
rotational instabilities may arise from non-radial toroidal 
modes e lmv (where m = ±1, ±2, . . . and ip is the azimuthal 
angle) . 

For sufficiently rapid rotation, the m — 2 bar mode be- 
comes either secularly or dynamically unstable. The onset 
of instability can typically be marked by a critical value of 
the dimensionless parameter /3 = T/|W|, where T is the 
rotational kinetic energy and W the gravitational potential 
energy. Uniformly rotating, incompressible stars in Newto- 
nian theory are secularly unstable to bar-mode formation 
when p > /3 SC c — 0.14. This instability can grow only in 
the presence of some dissipative mechanism, like viscosity 
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or gravitational radiation, and the associated growth time- 
scale is the dissipative time-scale, which is usually much 
longer than the dynamical time-scale of the system. By con- 
trast, a dynamical instability to bar formation sets in when 
P > Aiyn — 0.27. This instability is present independent 
of any dissipative mechanism, and the growth time is the 
hydrodynamic time-scale. 

In addition to the classical dynamical instability men- 
tioned above, there have been several studies indicating that 
a dynamical instability of the rotating stars occurs at low 
T/\W\ , which is far below the clas sical criterion of dynamica l 
instability in Newtonian gravity. iTohline fc Hachisul <ll99Cf) 
find in the self-gravitating ring and tori that an m = 2 
dynamical instability occurs around T/|W| ~ 0.16 in the 
lowest case in the cent rally condensed configurations. For 
rotat ing stellar models, IShibata. Karino fc Eriguchil <|2002| . 
2003) find that m — 2 dynamical instability occurs around 
T/\W\ ~ O(10~ 2 ) for a high degree (n c /D. s « 10) of differ- 
ential rotation. Note that Q c and fi s are the angular veloc- 
ity at the centr e and a t the equatorial surface, respectively. 
ICentrella et alJ (1200 ll) found dynamical m = 1 instability 
around T/\W\ ~ 0.09 in the N = 3.33 polytropic 'toroidal' 
star with a high degree (Q c /Q s = 26) of di fferential rota- 
tion, and lSaiio. Baumgarte fc Shapiro! <l2003t) extended the 
results of dynamical m = 1 instability to the cases with 
polytropic index iV > 2.5 and Q, c /Q s > 10. 

Computation of the onset of the dynamical in- 



© 2006 RAS 



2 M. Saijo and S. Yoshida 



stability, as well as the subsequent evolution of an 
unstable star, performed in a fully nonlinear hy- 
drodynamic simulation in N e wtonia n gravity, (e.g. 



Tohlinc. Durisen & McCollouehl 1 19851 


Durisen et alJ 


1983: 


Williams & Tohlinelll988l: iHouser. Centrella fc 


Smith 




ISmith. Houser fc Oentrellalll995 


; 

c 


Houser & Centrella 


199(| 


IPickett. Durisen & Davis! Il99d 


Toman et alJ 


1998; 


New. 


Centrella & Tohline 


2000) hav 


shown that 


/3dyn 



depends only very weakly on the stiffness of the equation 
of state. Once a bar has developed, the formation of a 
two-arm spiral plays an important role in redistributing 
the angular momentum and forming a core-halo structure. 
/3dyn are smal ler for stars with high degree of differential 
rotation jTohline fc Hadiisul Il99d: IPickett et all Il996t 
IShibata et_aL[ ]2002. 200 3|). Simulations in rela ti vistic grav- 
itatio n IShibata. Baumgarte fc Shapiro! l2000t ISaiio et alJ 
l200ll) have shown that /3d yn decreases with the compaction 
of the star, indicating that relativistic gravitation enhances 
the bar mode instability. 

Recently, several studies have been focused on the col- 
lapse of the rotating stars leading to non-axisymmetric dy- 
namical in stabilities in th ree-dimensional hydrodynamics. 
iDuez. Shapiro fc Yol 120041) investigated the collapse of a 
differentially rotating N = 1 polytropic star in full gen- 
eral relativity by depleting the pressure and found that the 
collapsing star for ms a torus which fragments into nonax- 
isymmetric clumps. IShibata fc Seki guchi (2005J) investigated 
rotational core collapse in full general relativity and found 
that a burst type of gravitational waves was emitted. In 
addition, they argued that a very limited window for the 
rotating star satisfies to exceed t he threshold of dynami- 
cal instability in the core collapse. ISaiio! d2005h studied the 
spheroidal and toroidal configuration of the collapsing star 
in conformally flat gravity, and found that toroidal config- 
uration has a potential source of gravitational waves due 
to the e nhancement of the non-axisymmetric dynamical in- 
stability. IZink e t al. ( 2005) presented a fragmentation of an 
N — 3 toroidal polytropic star to both one-armed spiral 
and a binary system in full general relativity, depending on 
the type of initial density perturbation. There the authors 
confirm that the instabilities found in Newtonian gravity 
also a ppear in general r elativistic stars of astrophysical rele- 
vance, [o^^^ij] i2005l) performed gravitational collapse of 
unstable iron cores at the centre of evolved massive stars 
in Newtonian gravity. Their simulations contained the evo- 
lutions from implosion of iron core (the computations done 
in two-dimensional code) up to the post bounce phase, in 
which they found growth of unstable m = 1, 2 oscillations. 

One of the remarkable features of these low T/|W| in- 
stabilities is an app earance of the corotation modes . As 
it is pointed out by IWatts. Andersson fc Jones! <l2005h the 
low T/\W\ unstable oscilla tion of bar-typed one found by 
IShibata et alJ J2002L l2003Tl has a corotation point. Here, 
corotation means the pattern speed of oscillation in the 
azimuthal direction coincides with a local rotational an- 
gular velocity of the star. It is well-known in the con- 
text of stellar or gaseous disk system that the corotation 
of oscillation may lead to instabilities. For instance, there 
have been several density wave models proposed to ex- 
plain spiral pattern in galaxies, in which wave amplifica- 
tion at the cor otation radius of spiral pattern is a key is- 
sue Jshu1ll992h . Another example of importance of coro- 



tation is found in the theory of thick disk (torus) around 
black holes. Initiated by a discovery of a dynamical insta- 
bility of geometri cally thick disk by iPapaloizou fc Pringld 
( 1984, 1985, 1987), several authors have studied these insta- 
bilities jBlaeslll985allbl lDrurvlll985t iBlaes fc Glatzel 1986t 
Goldreich. Goodman fc Naravanlll986l : iKojimal Il986t Il989l : 
lGlatzell987allbHGoodman fc Naravanll98a) . Instabilities of 
these systems are thought not to be unique in their origin 
and in their characteristics. Some seem to be related to local 
shear of flow and to share a nature with Kelvin-Helmholtz 
instability. Others may be related to corotation of oscilla- 
tion modes with averaged flow on which the oscillation is 
present. The mechanisms of instabilities by corotation, how- 
ever, seem not unique. A s is remini scent to 'Landau ampli- 
fication' of plasma wave JStbd ll992). a resonant interaction 
of corotating wave with the background flow (in the case 
of Landau amplification, background flow is that of charged 
particles) may amplify the wave, by direct pumping of en- 
ergy from background flow to the oscillation. The other may 
be an overreflection of waves at the corotation w hich may 
be se en in waves propagating towards shear layer dAchesonl 
Il976ft . 

The main purpose of this paper, in contrast to the pre- 
ceding studies of this issue, is to investigate the nature of low 
T/[W| dynamical instabilities, especially to study the qual- 
itative difference of them from the classical bar instability. 
As is mentioned above, recent studies have shown that dy- 
namical instabilities are possible for different region of the 
parameter space of rotating stars. Observing the existence of 
dynamical instabilities whose critical r/| VK| value are well 
below the classical criterion of bar instability, it is natural 
to raise a question on whether these two types, 'high T/\ W\' 
and 'low T/|W|', of dynamical instability are categorized in 
the same type of dynamical instability or not. 

Our study is done with both eigenmode analysis and hy- 
drodynamical analysis. A non-linear hydrodynamical simu- 
lation is indispensable for investigation of evolutionary pro- 
cess and final outcome of instability, such as bar formation 
and spiral structure formation. The nature of instability as 
a source of gravitational wave, which interests us most, is 
only accessible through non-linear hydrodynamical compu- 
tations. On the other hand, a linear eigenmode analysis is 
in general easier to approach the dynamical instability of 
a given equilibrium and it may be helpful to have physical 
insight on the mechanism and the origin of the instability. 
Therefore, a linear eigenmode analysis and a non-linear sim- 
ulation are complementary to each other and they both help 
us to understand the nature of dynamical instability. 

As a simplified system mimicking the physical na- 
ture of the differentially rotating fluid, we choose to study 
self-gravitating cylinder models. They have been used to 
study general dynamical nature of such gaseous masses 
as stars, accretion disks and of stellar system as galax- 
ies. Although there is no infinite-length cylinder in the 
real world, it is expected to share some q ualitative simi- 



larities with realistic astrophysical objects (Ostrikcr 


196E 


|Robelll979l lBalbinskilll985l llshibashi & Andd 


19851 


198€ 


IBlaes & Glatzel Il986l: lGlatzellll987allrl iLuvtenl 


1988. 


1989 


1990; Cook. Shapiro fc Stephans! '20031. Especially it has 



served as a useful model to investigate secular and dynami- 
cal instabilities of rotating masses. These works took advan- 
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tage of a simple configuration of a cylinder compared with 
a spheroid. 

This paper is organized as follows. In Sections and |^] 
we present our hydrodynamical and eigenmode analysis re- 
sults of dynamical one-armed spiral and dynamical bar in- 
stability. We present our diagnosis of dynamical m = 1 and 
m — 2 instabilities by using a canonical angular momen- 
tum in Section and summarize our findings in Section 
|S] Throughout this paper we use gravitational units with 
G = 1. Latin indices run over spatial coordinates. 



choose the same data sets as ISaiio et alJ i2003l) for inves- 
tigating low T/[W| dynamical instabilities in differentially 
rotating stars (m odels I and II in Table U corresponds to 
Tables II and I of lSaiio et al.l <l2003t) . respectively). We also 
construct an equilibrium star with weak differential rota- 
tion in high T/|W|, which excites the sta ndard dynamical 
bar instability, (e.g.. Ichandrasekharlll969l) . The character- 
istic parameters are summarized in Tabled 

To enhance any m = 1 or m — 2 instability, we dis- 
turb the initial equilibrium mass density p cq by a non- 
axisymmetric perturbation according to 



2 HYDRODYNAMIC SIMULATIONS IN 
DIFFERENTIALLY ROTATING STARS 

Here, we briefly describe the basic equation of the per- 
fect fluid hydrod ynamics in Newtonian gravity. We follow 
ISaiio et alJ i2003t) to perform our 3 dimensional Newtonian 
hydrodynamics assuming an adiabatic T-law equation of 
state 



P = (T - l)pe, 



(2.1) 



where P is the pressure, F the adiabatic index, p the mass 
density and e the specific internal energy density. For per- 
fect fluids, the Newtonian equations of hydrodynamics then 
consist of the continuity equation 



dp djpv 1 

at dx i 



= 0, 



the energy equation 
de d(ev J ) 1 

dt + ~~dxT~ = ~r e 



and the Euler equation 

d(pvi) d(pv t v 3 ) 



-(r-i) 



Pvi, 



dx 1 ' 



dt 



dxi 



d(P + P vi 
dx* 



dx 1 



(2.2) 



(2.3) 



(2.4) 



Here v 1 is the fluid velocity, $ is the gravitational potential 



A$ = 4np, 

and e is defined according to 
e = (p £ ) 1/r - 



(2.5) 
(2.6) 



We use the same type of artificial viscosity pressure P v i B in 
ISaiio et alJ i2003T) . When evolving the above equations we 
limit the stepsize At by an appropriately chosen dynamical 
time. 

We construct differentially rotating equilibrium stars 
based on lHachisul |l986j ). We assume a polytropic equation 
of state only to construct an equilibrium star as 



Kp 



1+1 /N 



(2.7) 



where k is a constant, N is the polytropic index. We also 
assume the 'j-constant' rotation law, which has an exact 
meaning in the limit of d — > 0, of the rotating stars 



Jo 



r _.- (2-8) 

where fi is the angular velocity, jo is a constant parame- 
ter with units of specific angular momentum, and zu is the 
cylindrical radius. The parameter d determines the length 
scale over which Q changes; uniform rotation is achieved in 
the limit d — > oo, with keeping the ratio jo/d 2 finite. We 



p = Poq 



1 + 5 ME+1 + (5 P)£!_LJ>! 



(2.9) 



where _R oq is the equatorial radius, with 6^ = w 1.7 — 
2.8 xlO" 3 in all our simulations. We also introduce 'dipole' D 
and 'quadrupole' Q diagnostics to monitor the development 
of m = 1 and m = 2 modes as 



D 



I imtp\ 

(e r ) 

\ It 

I imtp \ 

(e v ) 



1 

M 
1 

~M 



p^±^dV, 

ZD 



(2.10) 



(2.11) 



respectively. 

We compute approximate gravitational waveforms by 
using the quadrupole formula. In the radiation zone, 
gravitational waves can be described by a transverse- 
traceless, perturbed metric hj^ with respect to flat space- 



time. In the q uadrupole formu la, 
jMisner. Thorne fc WheeleJll973h 



1 1 



2(P_ 

rdV 13 



TT 



is found from 



(2-12) 



moment of the mass distribution, and where TT denotes 
the transverse-traceless projection. Choosing the direction 
of the wave propagation to be along the rotational axis (z- 
axis), the two polarization modes of gravitational waves can 
be determined from 



2 (hxx 



h> 



For observers along the rotation axis, we thus have 



rh+ 
~M 

rh x 
M 



(i tt 



2M dt 2 

1 d xTT 

M~d¥ xy 



i TT ) 

± yy J' 



(2-13) 

(2-14) 
(2-15) 



The time evolutions of the dipole diagnostic and the 
quadrupole diagnostic are plotted in Fig. We determine 
that the system is stable to m = 1 (m = 2) mode when 
the dipole (quadrupole) diagnostic remains small through- 
out the evolution, while the system is unstable when the 
diagnostic grows exponentially at the early stage of the evo- 
lution. It is clearly seen in Fig. Q that the star is more un- 
stable to the one-armed spiral mode for model I, and more 
unstable to the bar mode for models II and III. In fact, both 
diagnostics grow for model I. The dipole diagnostic, however, 
grows larger than the quadrupole diagnostic, indicating that 
the m = 1 mode is the dominant unstable mode. 
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Table 1. Three differentially rotating equilibrium stars that trigger dynamical instability. 



Model 


N a 


d/ Req b 




£~2 C 1 f2 s 


Pc/Pmax 


^maxd/^eq 


T/\W\ 9 


Dominant unstable mode 


I 


3.33 


0.20 


0.413 


26.0 


0.491 


0.198 


0.146 


m = 1 


11 


1.00 


0.20 


0.250 


26.0 


0.160 


0.383 


0.119 


m = 2 


111 


1.00 


1.00 


0.250 


2.0 


0.837 


0.388 


0.277 


m = 2 



a N: polytropic index. 6 i?eq: equatorial radius. °R p y. polar radius. d fl c : central angular velocity; Q s : equatorial surface 
angular velocity. e p c : central mass density; p max : maximum mass density. * R max d : radius of maximum density. 9 T: 
rotational kinetic energy; W: gravitational binding energy. 
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Figure 1. Diagnostics D and Q as a function of t/P c for three 
differentially rotating stars (see Table HI . Solid and dotted lines 
denote the values of D and Q, respectively. The Roman numeral in 
each panel corresponds to the model of the differentially rotating 
stars, respectively. Hereafter P c represents the central rotation 
period. 



The density contour of the differentially rotating stars 
are shown in Fig.|5]for the equatorial plane and in Fig.EJfor 
the meridional plane. It is clearly seen in Fig. H that one- 
armed spiral structure is formed at the intermediate stage of 



the evolution for model I, and that bar structure is formed 
for models II and III once the dynamical instability sets in. 

We show velocity fields in Fig.0]in the equatorial plane 
and in Fig. |S] in the meridional plane during the evolution. 
Note that shocks occur during the formation of m = 1 insta- 
bility. We also find that the fluid motion of the z-direction 
does not play a dominant role in generating the dynamical 
instabilities. 

We also show gravitational waves generated from dy- 
namical one-armed spiral and bar instabilities in Fig. |5] For 
m — 1 modes, the gravitational radiation is emitted not by 
the primary mode itself, but by the m = 2 secondary har- 
monic which is simultaneously excited, albeit at the lower 
amplitude. Unlike the case for bar-unstable stars, the grav- 
itational wave signal does not persist for many periods, but 
instead damp fairly rapidly. 



3 STABILITY ANALYSIS OF A 

DIFFERENTIALLY ROTATING CYLINDER 

3.1 Rotating selfgravitating cylinder model 

lOstrikeJ <ll965l) numerically computed physical characters of 
non-rotating infinite cylindrical masses. We here follow his 
treatment and introduce the normalization of variables. 

A cylinder rotates around its axis with a given angular 
velocity profile, which is a function of a cylindrical radial co- 
ordinate m. An equilibrium of the cylinder is determined by 
the balance between pressure gradient, centrifugal force and 
self-gravity of the cylinder. We also introduce an azimuthal 
angular coordinate ip, and a 2— coordinate set along the axis 
of the cylinder. 

For a fluid equation of state, we assume to have a poly- 
tropic relation 



P = Pc 



' = Pc 



3iV+l 



(3.1) 



where p c and p c are the normalization factors for a mass 
density and a pressure, which we choose those values on the 
rotational axis of the cylinder, N is a polytropic index. A 
variable with a bar denotes a normalized one. As shall be 
seen below, it is possible to construct a cylinder with a finite 
cylindrical radius by using large N. A non-rotating spheroid 
with N ^ 5 has an infinite radius, while a cylinder with 
TV = 25 still has a finite cylindrical radius. This qualitative 
difference from the well-known characteristics of spheroid 
simply results from the difference of geometry. 

We normalize the cylindrica l radial coordinate as zu — 
aw, where a = yf (N + 1)/AtyG, and the rotational angu- 
lar frequency as £7 = Q^/inGpc- Following a similar pro- 
cedure to obtain Lane-Emden equation of spherical poly- 
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Figure 2. Density contours in the equatorial plane for three 
differentially rotating stars (see Table [Tt . Models f, II, and III 
are plotted at the parameter (t/P c , Pmax/Pmax) = (16.2, 3.63), 
(134, 1.26), and (5.49, 1.20), where Pmax is the maximum den- 
sity, Pmax is the maximum density at t = 0, and R is the 
equatorial radius at t = 0. The contour lines denote densi- 
ties p/p max = 10-( 16 - l ) x0 ' 287 (j = 1,-..,15) for model I and 
p/p max = 6.67(16 -i) X 10~ 2 (i = 1, • ■ • , 15) for models II and III, 
respectively. 
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Figure 3. Density contours in the meridional plane for three 
differentially rotating stars (see Table . The parameters and 
the contour levels are the same as Fig. [5] 



tropes, more specifically, with using vj— component of the 
Euler equation and Poisson equation for a gravitational po- 
tential, we obtain Lane-Emden equation for a differentially 
rotating cylinder as 



d 1 



1 d6 



.-2 

dvj z vo avo 



2 nA 

dvj 



(3.2) 



The rotational profile of an angular velocity we study here is 
the same as the one used in the hydrodynamical simulations 
(equation 12. 811 . 



B 



+ A 1 



(3.3) 



where A and B are parameters. This is the same as equa- 
tion 12. 81 . with d = \f~A and jo = B. For simplicity, we 
hereafter omit 'bars' from all the equations. 

A frequently used dimensionless measure of rotation is 
r/| W| . The rotational kinetic energy T is defined as 



t = i / P zu 2 n 2 dv, 



(3.4) 



where integration is done for cylinder of unit le ngth. As to 
gravit ational energy, we follow the definition in lCook et alJ 
(2003) as 



w ■ 



pnV^dV = 



4-k 



(3.5) 



Here, $ is the gravitational potential and n ! is a unit normal 
vector of w =const. surface. The integration is performed 
for a unit length along the axis. m(ro) is the mass contained 
inside the cylindrical radius per unit length. 

For a large degree of differential rotation, the density 
maximum of the configuration becomes off-centred. An ex- 
ample of the profile of Lane-Emden function in such a case 
is plotted in Fig. [7| 



3.2 Linear perturbation of a cylinder 

3.2.1 Perturbation equations 

To study linearized oscillations of selfgravitating cylinder, 
we simultaneously solve linearized version of (1) equation 
of continuity; (2) Euler equation; (3) Poisson equation for a 
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Figure 4. Velocity field (v z /\v l B ^|) in the equatorial plane for 
three differentially rotating stars (see Table 0. The time for each 
snapshot is the same as in Fig. [2] Note that the velocity arrows 
are normalized as indicated in the upper right hand corner of each 
snapshot. \v\ | denotes the surface absolute velocity at t = 0. 
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Figure 5. Velocity field (v 1 /\v\ |) in the meridional plane for 
three differential rotating stars (see Table 111 . The time for each 
snapshot is the same as in Fig. |21 




Figure 6. Gravitational waveform for three differentially rotat- 
ing stars (see Table HI as seen by a distant observer located on 
the rotational axis of the equilibrium star. 



© 2006 RAS, MNRAS in press. [TUTU 



Low T/\W\ dynamical instability 7 



1.2 
1 

0.8 

<= 0.6 
0.4 
0.2 


0.2 0.4 0.6 0.8 1 

Figure 7. Equilibrium profile of Lane-Emden function for a 
polytropic index N = 25 and parameters of rotation (A, B) = 
(0.43, 1.6). Note that T/\W\ = 0.452. We find an m = 1 unstable 
mode in corotation in this case. 
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Figure 8. Imaginary part of the eigenfrequency a for a dynami- 
cally unstable m = 1 mode in a cylinder with A = 0.6. Note that 
a polytropic index N is 25. The parameters of models marked as 
(a) to (d) are listed in Table EJ 




gravitational potential. The adiabatic index of a perturbed 
fluid is assumed to coincide with that of background one, 
so that no p-mode appears in our computation. We also 
assume that there is no motion along the rotational axis of 
the cylinder, and therefore no dependence of the quantities 
on z-coordinate. 

Assuming a simple harmonic dependence of Eulerian 
perturbation of variable / as 

Sf(t,zu,ip) = Sfi (zu) exp(— iat + imip), (3-6) 

we can write down the perturbed equations as follows. 
Equation of continuity: 

du Ns I N d9 1 \ m» „ 
dzu u \ u dzu w J zu 
tu-component of Euler equation: 

su + 2 nv=-p- + -p-. (3.8) 
dzu dzu 

(^-component of Euler equation: 



(3.9) 



where k := ^2^(20, + zudfl/dzu) is the epicyclic frequency. 
Poisson equation for a gravitational potential 



jlH j oV = Ne 1- 

dzu z zu dzu zu 



(3.10) 



We have here defined Eulerian perturbation quantities as 
u = i • 6vY , v = 5vf, y — q = 59\. (3-H) 

We also define the following quantity 

s(zu) = a — mQ(za). (3-12) 

If s = for a certain cylindrical radius, the equations be- 
comes singular there. We call it corotation singularity. In 
that corotation radius zu CI t is defined as this singular- 

point. Corotation radius corresponds to a cylindrical sur- 
face on which the pattern speed of the oscillation 5R[<r/m] is 

1 Unless the rotational angular frequency is pathological, we ex- 
pect to have a regular singularity there. 



equal to the local angular frequency of the background flow, 
f2(tJ7 C rt)- Although s = is satisfied only for a purely real 
eigenfrequency, we here denote that a mode is in corotation 
when a cylindrical radius satisfies 3?[<r] — mQ, = in the 
cylinder, even if we have a complex eigenfrequency a. Thus 
a corotation radius is also defined for complex mode where 
the equation has no singularity. 

To close the eigenvalue problem, we should impose 
boundary conditions. One of the two conditions at the sur- 
face zu = zu s where equilibrium pressure becomes zero is the 
conventional free boundary condition. This means no stress 
is exerted on the cylindrical surface, which reduces to the 
condition: 

sq+^-u = Q. (3.13) 

dzu 

The other is the condition on the perturbed gravitational po- 
tential. From equation 13.101 . the perturbed potential out- 
side the cylinder (9 = 0, q = 0) is y ~ zu ±m . At the surface of 
the cylinder we impose the condition that the gravitational 
potential smoothly matches to the non-diverging solution at 
infinity. Therefore, the continuity condition of the potential 
requires 



dy m 

■j L + -y 



0, 



(3.14) 



at the cylindrical surface. On the rotational axis of the cylin- 
der, zu = 0, we impose a regularity condition on the vari- 
ables. 

Our numerical method to solve this eigenvalue problem 
is straightforward: We use the conventional shooting method 
to find eigenmodes. The system of linearized equations are 
solved both from the rotational axis and from the surface 
of the cylinder, with the boundary conditions being taken 
into account. At the intermediate matching radius, we im- 
pose a condition that both solutions connect smoothly. This 
procedure picks up the physical solutions of the eigenvalue 
problem. 

As we are interested in a dynamical instability, we 
assume the eigenfrequency takes complex values as well 
as other perturbed variables. We focus to the case where 
s(zu) = a — Q(zu) is non-zero, which is true except for a real 
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Figure 9. Eigcnfunction of an unstable m = 1 mode in N = 25 cylinder. We normalize the cylindrical radius zu with the surface radius 
vj s - Real part of Lane-Emden function q = 56 are plotted. Note that the normalization of the eigenfunction is arbitrarily. The label (a) 
to (d) in the figure correspond to the model in Table [5] 



a at the corotation region, i.e., Q.(w = zn s ) ^ o ^ SI(zj — 0). 
2 As a result, our present code can compute modes without 
corotation singularity. We also find convergence to 'modes' 
whose frequency is in the corotation region on the real axis 
of a plane. Although the 'eigenfunction' of it suggest that 
we are picking up one of the singular eigenmodes from con- 
tinuous spectrum, we cannot prove it by our present code. 

3.2.2 Characters of m — 1 unstable etgenmode 

First we note that we find unstable m = 1 mode only for ex- 
tremely soft equation of state. In general, we can construct 
a cylinder with a finite cylindrical radius (but an infinite 
length for z-direction) for an extremely large polytropic in- 
dex N than the case of a spheroid. 

From the r esults of hydrodynamical simulations 
JSaiio et aUl2003h . we find that an m = 1 instability ap- 
pears in a soft equation of state (N > 2.5). However in case 
of a cylinder, iV should satisfy iV > 20 to find an m = 1 
instability. This is a drawback of the cylinder to be com- 
pared with a result of the spheroid, although the behaviour 
of the unstable mode is similar to the spheroid. We expect 

2 The pure corotation needs careful treatment around the sin- 
gular point in order to pick up solutions as regular as possible 
We can use for instance Frobenius expansion <Ruoff fc Ko kkotas 
l200lHWatts et alj|2003f> there. 



that a qualitative nature of these modes are similar even the 
polytropic indices are quite different. 

For a fixed polytropic index, we have two parameters A 
and B to specify an equilibrium model. We construct a se- 
quence fixing A with changing B, which roughly corresponds 
to changing !T/|W| in a fixed degree of differential rotation. 
We summarize the characteristics of models studied in this 
paper in Table [5] 

In Fig. |H] we plot an imaginary part of the eigenfre- 
quency of m = 1 mode as a function of !T/| W|, fixing A = 0.6 
and iV = 25. The mode has a corotation radius inside the 
cylinder. The unstable mode appears only in a limited range 
of T/| W|. This behaviour of the imaginary part of the eigen- 
frequency is almost insensitive to the degree of differential 
rotation which is parametrized by y/A. 

Interestingly, this is reminiscent to the character of low- 
T/ \W\ bar (m = 2) i nstabi lity which is r ecently found 
bv IShibata et alJ (l200l l2003t) (see Fig. 3 in IShibata et ail 
(2003)). For a fixed degree of differential rotation, T/|W| 
that permits unstable mode is limited in a finite range. This 
may suggest that m — 1 instability and low-T/| W\ bar insta- 
bility may have the same (or similar) generation mechanism. 

We note that a real part of eigenfrequency is mono- 
tonically increasing function of T/|W| and of the degree of 
differential rotation (see Tabled . The frequency is the order 
of unity for all cases here. For each equilibrium model with 
a sufficiently large degree of differential rotation and a lim- 
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Table 2. Parameters for equilibrium models and eigenfrequency of m = 1 
and m = 2 modes. Models (a)-(d) correspond to those in Fig. l9l 



Model 


mode 


Af 




T/\W\ 






(a) 


m = 1 


25 


13.96 


0.4543 


0.5126 + 0.01726i 


0.2709 


(b) 


m = 1 


25 


11.34 


0.4597 


0.5512 + 0.03153i 


0.2806 


(c) 


m = 1 


25 


8.218 


0.4675 


0.6281 + 0.03982; 


0.2864 


(d) 


m = 1 


25 


6.067 


0.4722 


0.6910 + 0.01716i 


0.2715 


(b)-sl 


m = 1 


25 


11.34 


0.4597 


-0.245 (real) 




(b)-s2 


m = 1 


25 


11.34 


0.4597 


1.15 (real) 




(e) 


m = 2 


1 


13.00 


0.170 


0.3269 + 0.01256? 


0.507 



a a: eigenfrequency of the mode. b zu cr f. corotation radius; ra s : surface ra- 
dius. 



ited range of T/|VK|, we find unstable m = 1 modes in the 
corotation. Note that we also have exponentially damping 
modes, whose eigenfrequency are the complex conjugate of 
the dynamically unstable modes. We also find discrete sta- 
ble modes outside the corotation region, which have a purely 
real eigenfrequency. For equilibrium models where we found 
unstable m = 1 corotating modes, we did not find any dy- 
namically unstable m = 1 mode outside the corotation. 

We show the real part of the eigenfunction q of an m = 1 
unstable mode in Fig. [5] In order to compare the m = 1 
eigenfunction of the linear analysis with that of the hydro- 
dynamical simulation, we plot a perturbed mass density in 
a cylinder model (Fig. 1101 and a perturbed m = 1 unstable 
mass density in a differentially rotating star (Fig. lilt . In or- 
der to compute the perturbed m = 1 unstable mass density 
in a differentially rotating star, we follow 



Sp 



Pit) 



5p m=1 = - 



Peq, 

5p e~ 



°dip. 



(3.15) 
(3.16) 



We have a common feature that the m = 1 unstable mass 
density has a single oscillation inside the star and in the 
cylinder. However the behavior is quite different: for a cylin- 
der the oscillation is concentrated inside the corotation ra- 
dius, while for the differentially rotating star the oscillation 
is spread out to the whole equatorial surface radius. 

In order to focus on the comparison of the spiral struc- 
ture of the m = 1 unstable mode, we introduce a phase- 
constant curve of a perturbed radial velocity. In the linear 
analysis, complex velocity Sv^ is written as (we omit the 
factor of time dependence e _l<Tt since it is not relevant to a 
momentary spacial pattern) 

Sv m = U(vj)e lS(m)+zmv , (3.17) 

where U(w) is a real amplitude and S(w) is a phase func- 
tion. An equation defining phase-constant pattern is there- 
fore, 



S(vo) + mip — tan 



+ mip = const. 



(3.18) 



In a similar way, we also obtain a phase-constant pattern 
from the result of nonlinear hydrodynamical simulation. 
First we expand a perturbed radial velocity in the equatorial 
plane in terms of <p as 

i(mip-{-S r n) 



5v 



\\5v„ 



(3.19) 




Figure 10. Perturbed mass density of an unstable m = 1 mode 
for model (b) (see Table |5J. As it is a solution of linear problem, 
the scaling is arbitrarily chosen. 
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Figure 11. Amplitude of perturbed mass density of an unstable 
m = 1 mode for model I (See Table HI . Note that the amplitude 
is normalized with its maximum. Solid, dashed, and dash-dotted 
lines represent t = 8.15P C , 9.32P C , and 10.48P C , respectively. Ver- 
tical line represents the corotation radius of the rotating star. 



where 

8v™ = 



2^ ' SV ~ ' 



vj —imip 



dip, 



(3.20) 
(3.21) 
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Figure 12. Phase-constant curve of velocity perturbation 5v^ 
for model (b) (see Table |2J. The curve is plotted in the plane 
perpendicular to the rotational axis. The radius of cylinder is 
normalized to be unity. The dotted half-circle marks the corota- 
tion radius of the mode. Note that the direction of rotation of 
background flow is clockwise. 



S n 



= tan 



(3.22) 



Next, we focus on the phase-constant curve mip + S m = C, 
where C is a constant. Hereafter we choose C = since it 
only shifts the azimuthal angle. The m — 1 phase constant 
curve of the perturbed radial velocity in the equatorial plane 



VJ COS ifi 

vo sin (p - 



\Sv? 



\5v% 



■t*[6v?], 



(3.23) 
(3-24) 



We compare the phase-constant curve of m = 1 unstable 
mode between a cylinder (Fig. 1120 and a star (Fig. 1130. Both 
figures clearly have anm — 1 unstable spiral structure. How- 
ever the direction of arms are opposite. Also trailing or lead- 
ing nature of arms depends on the radial distance from the 
center in both cases. For the cylinder model in Fig. 1121 the 
arm changes from leading to trailing if we follow it from the 
center. For the star in Fig. 1131 an arm is leading inside while 
trailing outside. This apparent difference, however, does not 
prevent us from compa ring these t wo models. In fact, ac- 
cording to the result bv lRobd <ll979T) leading/trailing nature 
of spiral arm is not relevant to stability nature nor classifica- 
tion of eigenmodes. The same unstable eigenmode can have 
an arm with trailing, leading or mixed direction, depending 
on the equilibrium parameter. Thus the apparent difference 
is not significant here. The important point here is that the 
two model share one-armed spiral characteristics. 



4 CANONICAL ANGULAR MOMENTUM TO 
DIAGNOSE DYNAMICAL INSTABILITY 

4.1 Introduction of canonical angular momentum 

In order to diagnose the oscillations in a rotating fluid, we 
here introduce the canonical angular momentum following 



0.5- 



cr 



-0.5- 




-0.5 





x/R 



0.5 



eq 



Figure 13. Phase-constant curve of perturbed radial velocity in 
the equatorial plane of an unstable m = 1 mode for model I (see 
Table HI . The snapshot times are the same as that in Fig. 1111 
Solid circle denotes the corotation radius of the star. Note that 
the direction of rotation of background flow is counter-clockwise. 



iFriedman fc Schutzl <1978al) . In the theory of adiabatic lin- 
ear perturbations of a perfect fluid configuration with some 
symmetries, it is possible to introduce canonical conserved 
quantities associated with the symmetries. When we intro- 
duce Lagrangian displacement vector which is a vector 
connecting a perturbed fluid element to a corresponding 
non-perturbed one, the linearized perturbation equations 
are cast into a single second order differential equation for 
C (IFriedman fc SchutzllT978ah : 



a' ate + B)d t e + c)e = o. 



(4.1) 



Here the operators A and C are Hermitian, B anti- 
Hermitian, respecti vely up to divergence term [see equa- 
tions (32) to (34) in lFriedman fc Schutzl lll978al) for the pre- 
cise expressions], with respect to an inner product of dis- 
placements rf and £ l , 



(4.2) 



where rf means Hermite conjugate of r\. 

The master equation 14.111 is derived from the varia- 
tional principle with an action 



dt / dV £, 



where Lagrangian density C is denned by 



(4.3) 



(4.4) 



Applying Nother's theorem to t his Lagrang ian, we ob- 
tain canonical conserved quantities JWaldlll984l) . In partic- 
ular, if we have an axisymmetry in the background fluid, 
the corresponding Killing vector d v produces a conserved 
current related to the angular momentum, 

dC „ ^ dC 



J a = dZC 



-£ B e 



-£ P* 

a. . s j 



(4.5) 



where is a— component of vector 9^, and £ v denotes Lie 
derivative along a vector v. Note that our variable £ l is a 
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complex vector field. The time component of this current 
is a canonical angular momentum density, while spacelike 
component s are the flux density . In the rest of this section, 
let us follow lFriedman fc Schut3 (Il978al) to derive an explicit 
form of the canonical angular momentum. 

A natural symplectic structure is introduced as an inner 
product of a field £ and its conjugate momentum density, 



A{dtTj) + -Bri,Z). (4.6) 



It is easy to find that this product is conserved using the 
master equation 14. in and using the symmetric property of 
operators A, B, C as 



d t W(r,,0 = Q, 



(4.7) 



where both rf and are solutions to the master equa- 
tion jUJ . We obtain the canonical angular momentum of 
the system when the background fluid is axisymmetric (thus 
A,B,C commute with d v ), 



(4.8) 



which is a volume integral of the canonical angular momen- 
tum density defined by equation 14.511 . From the definition, 
these conserved quantities are quadratic in the Lagrangian 
displacement vector. 

These quantities, however, are 'gauge-dependent' in 
general. This means that we may find 'trivial displacement' 
vector £* for any physical solution of the master equa- 
tion 14.11 . The trivial displacement is added to 'physical' 
solution to produce a different displacement, which corre- 
sponds to the same physical solution (i.e., it does not change 
Eulerian perturbation of physical quantities). The trivial de- 
fines a class of gauge transformation of the same physical 
solution, under which the canonical energy or momentum 
are generally not invariant. Thus a naive use of it to the 
stabi lity problem of a fluid may lead to a wrong conclu- 
sion. iFriedman fc Schutzl lll978al) showed that we can find 
a class of a physical solution to the master equation 14.111 
called 'canonical displacement' orthogonal to the trivials, for 
which we have no contribution from the trivial displacement 
to the canonical quantities. Fortunately, in the case we are 
interested in (normal-mode problem with non-zero complex 
frequen cy), the displacement are a lways orthogonal to the 
trivials (Friedman fc Schut3ll978tJ) . 

Finally, we remark that the canonical energy or angu- 
lar momentum of the dynamically unstable modes are zero. 
This is simply found if we are conscious of the existence 
of an imaginary part aj in the eigenfrequency a. Then the 
conservation equation of the product W is written as 



d t W = 2aiW = 0, 



(4.9) 



and therefore W = 0. 

Next we write down the explicit form of canonical an- 
gular momentum used in the following discussion. 



Thus, 

MO = -\j d v £A)d t edV-\J d v i*B)d t edV 

(4.11) 

Note that since A, B, C are defined in the background fluid 
whose physical quantities are purely real, they are 'real' 
quantities (although B is anti-Hermite as an operator). As 
we are interested in a normal mode solution with harmonic 
dependence in t and <p, the displacement vector can be writ- 
ten as 



= Co Me 



in t — imip 



(4.12) 



The first and third terms in equation 14.111 are combined 
to produce 



(4.13) 



mfft[a] / p\£\ dV, 
Jv 

while the second and fourth are simplified as 



(4.14) 



> f P -Q[£v k V k C] dV. 
Jv 



Note that B is antisymmetric up to a divergence term which 
appears in the integral above. We have an exact cancellation 
to the contribution of B. As we are interested in the case 
with a circular flow as a background whose nonzero compo- 
nent of velocity is v v — Q(zu), we can easily find that 

6> fc Vkf = imfl\£\ 2 - zuQC'C + mSl^'C- (4-15) 

Finally, we get the simple form of the canonical angular 
momentum J c as 3 

J c = m I (JR[o-}-mQ.)p\i\ 2 dV-2rn I pzufl-Q[C '£"*]dV.(4.16) 
J v J v 



4.2 Application to oscillations of a cylinder 

We here present typical distribution of the canonical angular 
momentum in the cylinder model. The absolute amplitude 
of the plotted function here has no significance, since a linear 
eigenfunction can be scaled arbitrarily. 

In Fig. 1141 we plot the integrand of canonical angular 
momentum, defined from equation l|4.16^ as 

mj a (w) = m(K[cr] - mQ)p\i\ 2 - 2m,pvjQ, ■ ( 4 - 17 ) 

for unstable m — 1 mode. Here, j c is the canonical angu- 
lar momentum density. Note that an integral in the entire 
cylinder is zero for these cases. The features of the canonical 



-\(d^,A(d t + \Bt 
+\(Ad v dtt + \B^ 



(4.10) 



3 iNaravan et alJ il987Tl and lChristodoulou fc Naravar] Il992l) de- 
rived the same formula to study the oscillations of the slender an- 
nuli around a point mass gravity source. Watts (private commu- 
nication) derived a formula of the canonical energy of eigenmodes 
in differentially rotating shell of fluid. 
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Figure 14. Distribution of the canonical angular momentum density for m = 1 unstable mode. Plots are integrand of equation 14.161 . 
roj c (ro). Models (a) to (d) are the same as in Fig. [HI Vertical dashed lines mark the corotation radius of the mode. 



angular momentum distribution for m — 1 unstable modes 
are, 

(i) It changes sign around corotation radius zu CIt . 

(ii) It is positive for w < zu C rt, while negative for w > 
ro C rt- 

The feature |(i)| is remarkable and suggests us that the in- 
stability is related to the corotation. The feature has a clear 
contrast for a stable mode (Fig. 1151 . The canonical angu- 
lar momentum is either positive or negative definite, and it 
does not change its sign. Note that the former is the case 
when the pattern speed of mode is faster than the rotation 
of cylinder everywhere, while the latter is the opposite. This 
feature is expected from the equation (14.161 , if the first term 
is dominant. In such case, the sign of Sft[cr] — mQ(w) deter- 
mines the sign of the canonical angular momentum. This 
simple interpretation, however, does not hold for the dy- 
namically unstable mode. As it is shown in the feature f(ii)| 
above, we have a positive canonical angular momentum in- 
side the corotation, which is opposite to the sign of a — Cl(zu) 

for ^ W < TVcrt- 

In Fig. 1161 we show an example of canonical angular 
momentum distribution for m = 2 unstable mode of differ- 
entially rotating cylinder, w hich may be compared wit h the 
low T/\W\ bar instability of IShibata et alJ J200ll2005t . We 
did not find m = 2 unstable modes for the same param- 
eters as in the case of m — 1 instability. The features at 



13 




05/05. 



Figure 15. Distribution of the canonical angular momentum 
density for m = 1 stable modes. Plots are integrand of equa- 
tion 14.161 . roj c (ro). The parameters of top and bottom panels 
correspond to (b)-sl and (b)-s2 in Table 121 The model (b)-sl sat- 
isfies a — Q(zu) < 0, while the model (b)-s2 satisfies a — Q(zu) > 
throughout the fluid. 



the corotation radius, however, are the same as in m = 1 
instability. 

It is interesting to see how the profile of the canonical 
angular momentum changes when we consider the classical 
bar instability with uniform rotation. Unfortunately the bar 
mode of uniformly rotating cylinder has a neutral stability 



© 2006 RAS, MNRAS in press. [TUTU 



LowT/\W\ dynamical instability 13 



2 [ 

u 

B 

-2 - 



0.2 



0.4 



0.6 



0.8 



ro/tn 



Figure 16. Distribution of canonical angular momentum den- 
sity for m = 2 unstable modes (model (e) in Table |5J- Plots are 
integrand of equation 14.161 . mj c (m). The vertical dashed line 
marks the corotation radius. 



point at the breakup limit l|Luvtenl 1990). We instead looked 
at m = 2 instability o f unif ormly rotating, incompressible 
Bardeen disk <lBardeenlll975l) and the classical bar instabil- 
ity of Maclaurin spheroid (see Appendix[X]for these compu- 
tation). These are actually more suitable for comparison to 
differentially rotating spheroidal model, which we present 
in the following section. For both of the models we have 
analyt ic expressions of oscillatio n mode s [Schutz fcBo wenl 
( 1983) for a Bardeen disk and Chan drasekharl jl969T) for 
Maclaurin spheroid]. It is remarkable that the canonical an- 
gular momentum density is zero everywhere (which ensures 
that the total canonical angular momentum vanishes). This 
is in a clear contrast with the m = 2 instability in the cylin- 
der with highly differential rotation. 



4.3 Differentially rotating star 

We here present the application of the canonical angular 
momentum to our hydrodynamics results. Following three 
assumptions are made to adopt our hydrodynamic results 
of dynamical instability (for both m = 1 and m = 2) to the 
perturbative approach; 

(i) All physical quantities is in coherent oscillation of 
growing mode. 

(ii) For each model, a growing mode with single m is dom- 
inant. 

(iii) The motion of z-direction does not contribute to the 
instability. 

From assumption | (i) | we introduce a complex frequency 
o that represents the same growing mode m. Therefore all 
physical quantities f(t) that have a time dependence should 
satisfy 



/(*) = h exp(-io-t), 



(4.18) 



where a = or + iai , fi is a complex quantity. Note that or 
and <tt are real quantities. 

The assumption |(ii) | comes from the fact that single m 
mode has a dominant contribution to the dynamical instabil- 
ity in the diagnostics (Fig.0. Therefore all physical quantity 



Table 3. Eigenfrequency and the corotation radius of three dif- 
ferentially rotating stars 



Model 


(7 


[fie] 


ro cr t [-Req] 


I 


0.590- 


f 0.0896i 


0.167 


11 


0.284 - 


t- 0.0121i 


0.492 


111 


0.757 


+ 0.200i 





f(zu,p>) that have a dependence of azimuthal angle should 
satisfy 



/(•B7, (p) = f{vu) exp(imp) 
where 

1 



/(«0 = 



2tt 



dip f(w, ip) exp(-imip). 



(4.19) 



(4.20) 



From the velocity snapshots in the meridional plane 
(Fig. 0, the motion of the fluid along the rotation axis 
does not have a significant contribution to the instability. 
Therefore we can neglect the z-dependence in the function 
(assumption |(iii)} as 

f(m,z) = f(m). (4.21) 

From the three assumptions, we determine the fre- 
quency from the dipole and quadrupole diagnostics as 



D(t) 



1 

M 



dvp(t, w, ip, z) exp(i<p) 



exp(— iot)— J dzu2nzu p(zu), 



Q{t) = m I dvp(t,uj,tp,z)exp(2ip) 



exp(— iot)— J dvj2irwp(w). 



(4.22) 



(4.23) 



With the formulae above, we extract the complex fre- 
quency by fitting evolutions of dominant diagnostics for each 
cases. 

Note that we averaged the frequency in the early stage 
of the evolution, since the real part of the frequency is almost 
the same. Note also that D(0) ^ and Q(0) 7^ since 
we put m — 1 and m = 2 density perturbation at t — 
(equation 12.91 to initiate m — 1 or m = 2 dynamical 
instability. 

We summarize the corotation radius and the complex 
frequency in Table[jJ]from Fig.0 Eulerian perturbed velocity 
is defined as 



dtv % {t,w,p) = v l (t,m,tp) — Ve q (TO-), 



(4.24) 



where v z (t, tu, ip) is the velocity at t, vl^vj) is the velocity 
at equilibrium, respectively. Following assumption |(ii)| we 
find, 



dtv l (t, zu) 



1 

2^ 



dp d t v'(t, zu, p) exp(-imtp). (4.25) 



The Lagrangian displacement £ l should satisfy the fol- 
lowing differential equation: 



dtt, =o t v + £ a k v cq - v cq d k t; 



(4.26) 
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Using the three assumptions, vj and tp component of the 
Lagrangian displacement should be written as 

idtv™ 



r = 



e = 



— mf2 e 
id t v v 

— mQ c 



(a — mfi oq ) 2 



(4.27) 
(4.28) 



where f2 oq denotes the angular velocity at equilibrium. 

Using the Lagrangian displacement extracted by the 
formula above, we compute the canonical angular momen- 
tum density defined as an integrand of the canonical angular 
momentum in equation 14.161 as 



Jc 



where 



(4.29) 



wj c (w) = m(K[cr] - mQ) P \^\ 2 - ImpwO. ■ £*""]• ( 4 -30) 



Note that we have used the assumptions |(i) | - |(iii) | to compute 
the canonical angular momentum density. 

We show the snapshots of canonical angular momen- 
tum density in Fig. 1171 Since we determine the corotation 
radius using the extracted eigenfrequency and the angular 
velocity profile at equilibrium, the radius does not change 
throughout the evolution. For low T/|W| dynamical insta- 
bility, the distribution of the canonical angular momentum 
drastically changes its sign around the corotation radius, 
and the maximum amount of canonical angular momentum 
density increases at the early stage of evolution. This means 
that the angular momentum flows inside the corotation ra- 
dius in the evolution. On the other hand, for high T/|W| 
dynamical instability (the bottom panel of Fig. 1170 , which 
may be regarded as a classical bar instability modified by 
differential rotation, the distribution of the canonical angu- 
lar momentum is smooth and with no particular feature. 

Note that the amplitude of wj c is orders of magnitude 
smaller than those in the corotating cases in top and middle 
panels of Fig. 1171 Contrary to the linear perturbation analy- 
sis in Section ^. 21 the amplitude here is not scale free and the 
relative amplitude has a physical meaning. Thus the small- 
ness of it for model III suggest that it should be exactly zero 
everywhere in the limit of linearized oscillation. The devia- 
tion from zero may come from the imperfect assumption of 
linearized oscillation, that is made here to extract oscillation 
frequency and Lagrangian displacement vector. 

From these different behaviours of the distribution of 
the canonical angular momentum, we find that the mecha- 
nisms working in the low T/\W\ instabilities and the classi- 
cal bar instability may be quite different, i.e., in the former 
the corotation resonance of modes are essential, while the 
instability is global in the latter case. 



5 SUMMARY AND DISCUSSION 

We have studied the nature of three different types of dy- 
namical instability in differentially rotating stars both in 
linear eigenmode analysis and in hydrodynamic simulation 
using canonical angular momentum distribution. 

We have found that the low T/\W\ dynamical instabil- 
ity occurs around the corotation radius of the star by investi- 
gating the distribution of the canonical angular momentum. 
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Figure 17. Snapshots of the canonical angular momentum dis- 
tribution tuj c {tu) in the equatorial plane for three differentially 
rotating stars (see Table 0. Solid, dotted, dashed, and dash- 
dotted line represents the time t/P c =(3.47, 6.93, 10.40, 13.86) 
for model I, t/P c =(45.68, 56.43, 67.18, 77.97) for model II, and 
t/P c =(1.10, 2.19, 3.29, 4.39) for model III, respectively. Note 
that vertical line in panels I and II denotes the corotation radius 
of the star (model III does not have a corotation radius). Wc also 
enlarged the figure around the corotation radius for panels I and 
II, which is presented in the right upper part of each panel. Al- 
though our method of determining the corotation radius is not 
precise, we clearly find that the distribution significantly changes 
around the corotation radius. 



We have also found by investigating the canonical angular 
momentum that the instability grows through the inflow of 
the angular momentum inside the corotation radius. The 
feature also holds for the dynamical bar instability in low 
T/IW 7 ), which is in clear contrast to that of classical dynam- 
ical bar instability in high T/|W^|. Therefore the existence 
of corotation point inside the star plays a significant role of 
exciting one-armed spiral mode and bar mode dynamically 
in low T/|W|. However, we made our statement from the 
behaviour of the canonical angular momentum, the state- 
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ment holds only in a sense of necessary condition. In order 
to understand the physical mechanism of the low T/[W] dy- 
namical instability, we need another tool and it will be the 
next step of this study. 

The feature of gravitational waves generated from these 
instabilities are also compared. Quasi-periodic gravitational 
waves emitted by stars with m = 1 instabilities have smaller 
amplitudes than those emitted by stars unstable to the 
m = 2 bar mode. For m = 1 modes, the gravitational ra- 
diation is emitted not directly by the primary mode itself, 
but by the m = 2 secondary harmonic which is simultane- 
ously excited. Possibly this m = 2 oscillation is generated 
through a quadratic nonlinear selfcoupling of m = 1 eigen- 
mode . Remarkab l y the precedent studies dCentrella et alJ 
l200ll : ISaiio et alJ I2003T) found that the pattern speed of 
m = 2 mode is almost the same as that of m = 1 mode, 
which suggest the former is the harmonic of the latter. Un- 
like the case for bar-unstable stars, the gravitational wave 
signal does not persist for many periods, but instead is 
damped fairly rapidly. We have not understood this remark- 
able difference between m = 1 and m = 2 unstable cases. 
One of the possibility may be that the unstable m = 1 eigen- 
mode tends to couple to higher and higher m modes (which 
are not necessarily unstable and could be in the continuous 
spectrum) and pump its energy to them in a cascade way. 
However, we have not found the feature that prevents m = 2 
mode from this cascade dissipation. 

Another possibility is that the spiral pattern formed in 
m = 1 instability redistributes the angular momentum of the 
original unstable flow, so that the flow is quickly stabilized. 
Inside the corotation radius, the background flow is faster 
than the pattern, while it is slower outside. A similar mech- 
anism to Landau damping in plasma wave which transfer 
the momentum of wave to background flow may work at the 
spiral pattern. The pattern may decelerate the background 
flow inside the corotation and accelerate it outside the coro- 
tation, which may change the unstable flow profile to stable 
one. As we do not see a spiral pattern forming in the low 
T/|W| bar instability, it may eludes this damping process. 
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APPENDIX A: CANONICAL ANGULAR 
MOMENTUM OF CLASSICAL BAR 
INSTABILITY 

In this appendix we evaluate the angular momentum density 
of bar instability in Bardeen disk and in Maclaurin spheroid. 
This quantity is defined in the cylindrical coordinate as 

mK[s]p|£| 2 - 2mQp • fH 
P£ ($t{s]\8v\ 2 +2ftQ [8v m ■ Sv**]) , 



(Al) 
(A2) 



the general perturbation can be expanded by Legendre poly- 
nomials as 



oo oo 



2 , a i Pn+ m yn)e 



(A4) 



l — m— — oo 



Master equatio ns [equations (3.1) and (3.4) in 

ISchutz fc Bower] lll98aft ) for perturbation can be de- 
composed into those for each I and m. The bar mode 
is I = 0, m = c ase. F rom equations (3.13) to (3.15) 
of ISchutz fc Bower] dl983f) . we get a simple eigenvalue 
equation 



2 7T 

where Q is the angular velocity of disk 



8// 



(A5) 



(A6) 



while s is a mode frequency observed in a corotating frame 
with disk. Both of them are normalized by fl c \d, an angular- 
frequency of cold dis k, u is called 'aspect ra tio' parame- 
ter (equation (2.18) in lSchutz fc Boweni Jl983l) '). When this 
parameter is smaller than 7r/16, one of the modes above 
becomes dynamically unstable. The corresponding radial 
eigenfunctions are written as 



2 
«2 



if 1 + 



8F 



where 



r = 



s(s + 2Q.) 

sr 

s(s + 2D.) 



Phi), 



2/x 



(A7) 
(A8) 

(A9) 



16 7T 

Using a\ and /3|,Eulerian perturbation of velocity for bar 
mode is written as 



Sv^ — —i 



( — a\ + —02] = —QiKzu, 
V azu zu / 



d n2 2 2 CT s 

—— p 2 + -«2 = oKw, 
azu r 



(A10) 
(All) 



where we used the definition of r\. With bearing in mind that 
the real part of s is f2, equation 1A2I gives that j c is 0. 



where s — a — mfi is the frequency of the mode observed in 
corotating frame of the fluid. The components of Lagrangian 
displacement £ and Eulerian velocity 8v are written in the 
local orthonormal frame. 



Al Bardeen disk 

iBardeer] |[l975) analytically constructed a self-gravitating 
'warm' disk with finite thickness by finding corrections to 
an infinitesimally thin disk from which he started his ap- 
proximation. Per t urbat ion of a Bardeen disk is studied in 
ISchutz fc Bower] l)l983Tl . Here, we follow their definitions 
and notations. Velocity perturbation is written with two po- 
tential functions a and (5 as 



8v 



-iVa - V x (f3e z ), 



(A3) 



w here e z is th e unit vector in z-direction. Introducing r\ = 
— (vj/R) 2 in terms of cylindrical radial coordinate zu, 



A2 Maclaurin spheroid 

For the bar mode of Maclaurin spheroid, we have an analyt- 
ical expression of Lagrangian displacement, 

C = kzue ±tip , C = kzue ±l{Lp+ %\ £ z = 0, (A12) 

in Cartesian components llChandra sckhar 1969). Here fc is a 
constant. It is straightforward to see that the corresponding 
components in cylindrical coordinate are 



e 



kzue 



±2iip 



±ke 



±2iip 



0. 



(A13) 



£ v is defined with respect to coordinate base vector d/dtp. 
From equation <A1> with the fact that real part of frequency 
of unstable bar mode is il, we easily see that j c vanishes 
everywhere. 
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